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Transconductance and effective Lande factors for quantum point contacts: 
spin-orbit coupling and interaction effects 


K. Kolasihski, A. Mrenca-Kolasinska, and B. Szafran 
AGH University of Science and Technology, 

Faculty of Physics and Applied Computer Science, 
al. Mickiewicza 30, 30-059 Krakow, Poland 

We analyze the effective g* factors and their dependence on the orientation of the external mag¬ 
netic field for a quantum point contact defined in the two-dimensional electron gas. The paper 
simulates the experimental procedure for evaluation of the effective Lande factors from the transcon¬ 
ductance of a biased device in external magnetic field. The contributions of the orbital effects of 
the magnetic field, the electron-electron interaction and spin-orbit (SO) coupling are studied. The 
anisotropy of the g* factors for the in-plane magnetic field orientation, which seems counterintuitive 
from the perspective of the effective SO magnetic field, is explained in an analytical model of the 
constriction as due to the SO-induced subband mixing. The asymmetry of the transconductance as 
a function of the gate voltage is obtained in agreement with the experimental data and the results 
are explained as due to the depletion of the electron gas within the quantum point contact constric¬ 
tion and the related reduction of the screening as described within the DFT approach. The results 
for transconductance and the g * factors are in a good agreement with the experimental data [Phys. 
Rev. B 81 , 041303, 2010]. 


I. INTRODUCTION 


Quantum point contacts (QPC) - building blocks 
for quantum transport devices - turned out early 
to be spin-active elements when an enhancement of 
the effective Lande factor g* in QPCs for the two- 
dimensional electron gas confined to a quasi one¬ 
dimensional constriction was found (T] . The interac¬ 
tion effects [2], namely the exchange bias PHE] are 
held responsible for this enhancement as well as for 
the 0.7 conductance anomaly mm- In systems with 
strong spin-orbit coupling the QPCs act as spin-filters 
working in the absence of the external magnetic field 
CM7]. In the spin control and manipulation devices 
[IMS] the external magnetic field can be replaced 
by an effective momentum-dependent magnetic field 
m introduced by the spin-orbit (SO) interaction for 
moving carriers. The SO interaction results from the 
asymmetry of the structure or the crystal lattice and 
in turn introduces a spatial anisotropy of the spin phe¬ 
nomena, introducing in particular the conductance de¬ 
pendence on the orientation of the in-plane magnetic 
field |2()[I22| . The anisotropy of the g* factors for an 
in-plane and an out-of-plane orientation of the mag¬ 
netic field is clearly present in the experimental data 
[25] [2T . Additionally, the data of Ref. [25] exhibited a 
weak anisotropy of the g* factor for the in-plane ori¬ 
entation of the magnetic field i.e. perpendicular and 
parallel to the direction of the current flow. A stronger 


anisotropy of the g* factors for the in-plane orienta¬ 
tion of the magnetic field was observed for holes in 
quantum wires [2J>j and quantum point contacts [25 
and explained [271128] as due to the spin-orbit inter¬ 
action for the valence band carriers. The effective g* 
factors for electrons in quantum wires was studied in 
a number of papers |29l)52j . The combined effects of 
SO and electron-electron interactions for 0.7 anomaly 
in quantum point contacts were recently discussed [22 
for neglected mixing of subbands within the QPC con¬ 
striction. 

The purpose of the present paper is an analysis 
of the interactions and spin-orbit coupling effects for 
the effective g* factors in quantum point contacts de¬ 
fined within a two-dimensional electron gas. While 
the previous theories dealt with the Zeeman shifts of 
subbands PHI [251 EZM], we report on a numerical 
model focused on simulation of the standard experi¬ 
mental procedure ]25. (21} [551155] for extraction of the 
g* factors from the experimental transconductance. 
The procedure [23J EH EM6] exploits the transcon¬ 
ductance lines at the edges of half-quantized conduc¬ 
tance plateaux that appear in external magnetic field 
in presence of the source-drain bias. We analyze the 
anisotropy of the g* factor and its dependence on the 
subband index for Rashba SO interaction due to an 
external electric field. We find an unexpected reduc¬ 
tion of the effective g* factor for the external mag¬ 
netic field oriented in-plane and perpendicular to the 
current flow direction i.e. parallel to the nominal ori- 
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entation of the effective magnetic field due to the SO 
coupling. The procedure used for extraction of the g* 
factors from transconductance gathers the data of the 
bottom of subsequent subbands entering the transport 
window, for which the effective magnetic field is neg¬ 
ligible. We still find that the extracted g* factors are 
anisotropic for the in-plane orientation of the external 
field, and explain this feature by the SO coupling of 
the constriction subbands. 

We discuss the g* factor anisotropy with and with¬ 
out the electron-electron interaction. For the inter¬ 
action effects described within the spin density func¬ 
tional theory (SDFT) approach j5H5] we find a good 
agreement with the experimental g* factors of Ref. 
[25 ] . Moreover, we discuss the interaction effects due 
to removal of the electron charge density from the 
QPC constriction m and demonstrate that it leads 
to the asymmetry of the transconductance lines as a 
function of the gate voltage as observed experimen¬ 
tally [25;. We separate the contribution of the spin- 
orbit coupling effects and the orbital effects in the 
out-of-plane magnetic field. The electron-electron in¬ 
teraction alone induces a strong imbalance between 
the g* factor for the in-plane and out-of plane mag¬ 
netic field orientation. 




Figure 1. (a) Sketch of the considered system. A channel 
of width 240nm with QPC potential modeled by Eq. @. 
The channel is assumed infinite in the x direction with 
the computational box that covers 800 nm. The red line 
in the center of the QPC indicates the cross section of 
the system for which the dispersion relations in Fig. [2] 
were obtained, (b) Cross section of the system along the 
x-axis showing the potential distribution for a finite Vsd 
potential. We assume the symmetric drop of the potential 
along the device. The dashed line schematically shows the 
total potential profile Vsd + I4 x t- 


II. QPC MODEL AND TRANSPORT 
THEORY 

We consider a system defined within a two- 
dimensional electron gas of a wide channel contain¬ 
ing a constriction of a quantum point contact (QPC) 
in its center (see Fig. |lja)) . The channel’s width is 
240nm and its length is assumed infinite (computa¬ 
tional box covers 800 nm). The electrons flow from 
the input/source lead to the output/drain lead. We 
solve the scattering problem of the Fermi level elec¬ 
trons using the one-electron spin dependent effective 
band mass Schrodinger equation for Fermi energy Ep 


H is the momentum operator and A is 

the vector potential. The external potential I4 x t de¬ 
scribes the QPC narrowing and is assumed to be of 
Gaussian shape [35] modeled with following formula 

Rext = Vgate(x, y ; 400nm, 0, 50nm, 42nm) 

+ Vgat e (x, y; 400nm, 240nm, 50nm, 42nm) (2) 


where the upper and lower Vg ate terms correspond to 
the upper and lower QPC potential profiles (see Fig. 
[lja)) and are given by 


legate(^b y'l ®gi 2/gAx, A y) — Vge 


= We () e 


y-yg \ 2 

2Ay ) ' 


p 




eW, 


,3A*b B 


(1) 


+H lat + H 


rashba J 


= Ep 


where the plane of confinement is taken at z = 0 
and B = ( B x ,B y: B z ) is the magnetic field vector. 
We use the material parameters for the Ino. 5 Gao. 5 As 
alloy [31], the effective mass m e ff = 0.0465ttto and 
the Lande factor g = 9. In Eq. (1) p = hk = 


In the equation above (x g ,j/ g ) are the coordinates of 
the center of the Gaussian function with (A x , A y ) 
defining its width and height. The QPC constriction 
is controlled by the potential V g parameter which is 
an analogue to the electrostatic potential energy intro¬ 
duced by the gate [25]. The third term in the Hamilto¬ 
nian Q represents the Zeeman splitting for magnetic 
field B with <r being the Pauli matrices vector 
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We use the gauge A = ( eB y z — eB z y, — eB x z , 0). For 
the plane of confinement z = 0 the applied vector po¬ 
tential reduces to the one given by the Landau gauge 
A = (— eB z y 1 0, 0). Therefore, the in-plane compo¬ 
nents of the magnetic field B x , B y enter the Zeeman 
term only. On the other hand B z enters the mo¬ 
mentum operator and introduces orbital effects of the 
magnetic field to the kinetic energy and the SO inter¬ 
action. The latter is described by the last two terms 
of the Hamiltonian JT5] with the Rashba term 

-f^rashba = 'frsb x^y & y^x} (3) 

due to the electrostatic confinement of the 2DEG in 
the growth direction. The Rashba parameter 7 rs b = 
o^d-Fz defines the strength of the SO interaction in 
the given material, where « 3 D is a material constant 
and F z is the electric field in z direction. We set 
a 3 D = 0.572nm 2 as an average value for GaAs and 
InAs f40j. The values of the electric field in the growth 
direction are taken up to F z = 20 gO]. The value 
of F z can be controlled in the experiment by the gate 
voltage. We will discuss below the role of the 7 rs b 
parameter on the anisotropy of the effective g-factor. 
ifiat in Hamiltonian (1) describes the so-called lateral 
SO coupling present due to in-plane electric fields m 
with 


FLlat 



C>Kxt 

dy 


- fcv 


C>Kxt ] 

dx J 


This term is included for consistency, however its in¬ 
fluence on the effective g* factor turns out to be negli¬ 
gible due to i) the symmetry of the considered poten¬ 
tial with respect to the axis of the system m , ii) and 
small values of the in-plane gradients generated 
by the QPC potential. 

All the SO effects discussed in this paper result from 
the electric fields and we neglect the Dresselhaus SO 
interaction due to the asymmetry of the crystal lattice. 
The neglect of the crystal-asymmetry-induced SO in¬ 
teraction is justified for strong vertical electric fields 
- assumed below and for large width of the quantum 
well confining the electrons in the growth direction 

013 - 

In order to solve the scattering problem we use the 
finite difference implementation of the quantum trans¬ 
mitting boundary method (QTBM) [22, H3j for the 
Schrodinger equation 0 - A more detailed description 
of the matrix formulation of the problem and bound¬ 
ary conditions can be found in the Appendix A. Once 


the scattering problem is solved we calculate the con¬ 
ductance using the Landauer approach. In the linear 
response regime and for small temperatures T —> 0 
the formula can be put in the form 

p 2 p 2 M 

G=-T tot (E F )^-Y, T ^ ( 4 ) 

k 

where T\ is the transmission probability of the /c-th 
mode incoming from the input lead to the system, 
with M being the total number of transverse spinor- 
modes. In case of finite source-drain potentials we 
calculated the current with the formula 

+eVsD/2 

/(U SD ;T = 0 ) = ^ J T tat {E F + E)dE, (5) 

— eVsD/2 


where we assume a linear potential drop for non-equal 
chemical potentials of the leads. The sketch of the 
potential profile in the device when the finite Vsd po¬ 
tential is applied is depicted in Fig. 0 b)- The con¬ 
ductance is calculated by 


a/(Vs D ) 

<9EsD 


( 6 ) 


and the transconductance is next defined by dG/dV g . 

The finite difference discretization of the Hamilto¬ 
nian ([I]) was performed with grid spacing Ax = Ay = 
2nm and the Fermi energy E F = 50meV. For « 3 D = 0 
and \B\ = 0 the number of transverse modes in the 
input and output leads at the Fermi level was M =36. 


III. RESULTS WITHOUT THE 
ELECTRON-ELECTRON INTERACTION 

A. Dispersion relation 

In order to set the ground for the further anal¬ 
ysis we display in Figure [2] the dispersion relation 
and the spin polarization of separate subbands for a 
homogenous channel of lateral confinement that we 
adopted from the narrowest cross section of the QPC 
(red line in Fig. III. In the absence of the exter¬ 
nal magnetic held (Fig. [2ja-c)) the electron spins 
are polarized by the effective magnetic held result¬ 
ing from the Rashba interaction. In a general case 
the held is of the form H e // = 'yrsbiky, — k x , 0 ) (see 
Eq. 0)- In a channel with potential independent of 
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x, the wave vector k x = k is a well defined quantum 
number. On the other hand, k y is an operator that 
produces a vanishing contribution in terms of the first- 
order perturbation, due to a definite symmetry of the 
wave functions with respect to the axis of the channel. 
Hence, for a symmetric homogenous channel we have 
B e ff = 7 rs b(0, —k, 0). The dispersion relation of Fig. 
[2jb) indicates the Zeeman-like splitting of subbands 
with the lower-energy bands spin-polarized parallel to 
B e ff , i.e in the —y direction for k > 0, and +y direc¬ 
tion for k < 0. The SO coupling for B = 0 introduces 
shifts of the dispersion relations for y-polarized spin 
subbands along the wave vector scale [2], which are 
next shifted on the energy scale by the magnetic field 
oriented in the y direction, resulting in the asymme¬ 
try [20j of the dispersion relation of Fig. [2jg-i) . The 
polarization of subbands remains unchanged when the 
external magnetic field is applied parallel to the y axes 
[Fig. [2j[g-i)]. On the other hand for the magnetic field 
parallel to the x and z axis we observe a competition 
between the Zeeman effects due to the external field 
and the Rashba effective field. The latter prevails at 
high k. 


B. Transconductance and the g* factors 


We performed calculations to simulate the clas¬ 
sical procedure [23J 1551155] of extracting the 

g* factors from the experimental transconductance 
cPl/dVgdVsD data based on compensation of the Zee- 
man splitting by the source-drain bias. Figure [3j[a) 
shows the current as calculated from the Landauer 
formula ([5| for B = 0. The resulting conductance 
G = is displayed in Figure 3lb) and exhibits char¬ 
acteristic diamond shaped plateaux. The transcon¬ 
ductance of Figure [3^c) is non-zero only near voltages 
for which the number of subbands of the constric¬ 
tion that fall between the electrochemical potentials 
of source and drain changes. The plateaux of odd 
conductance (in units of Gq = %-') appear (35) when 
the number of conducting subbands for the +x and 
—x directions differ by one. A non-zero external mag¬ 
netic field [cf. Fig. [3](d)] introduces additionally Zee- 
man splitting of subbands producing half-plateux of 
conductance. 

The transconductance as a function of the voltages 
is used in the experiments HH US [33] [35] to extract 
the gate voltage to energy conversion factor. We use 
the data for B = 0 of Fig. [3jc) to evaluate the con- 
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Figure 2. Dispersion relation for a homogenous chan¬ 
nel of the lateral confinement taken from the QPC con¬ 
striction (red line in Fig. [l]). Results were obtained 
for q; 3 d = 0.572nm 2 and F z = 20 meV/nm which gives 
7 rsb = 11.4 meVnm. The columns indicate the average spin 
polarization in x , y and z directions (green color stands 
for (si) — 0, the red for (s;) = h/2, and the blue for 
(s») = —h/2). The first row (a-c) is for B = 0. The other 
plots correspond to B\ = 1.0 T, oriented in x direction 
(d-f), y (g-i) and z (j-1). 


version factor 

_ 1 rfVsp 
conv “ 2 dV g ' 

The 1/2 factor accounts for the shift of the source 
(drain) potentials, which is half of the applied bias: 
l/2eV sd (—l/2eV s d). The derivative is evalu¬ 
ated as the ratio of the source drain voltage to the 
A V g splitting of the transconductance lines for the 
lowest subband [see Fig. Sc)]- The conversion fac¬ 
tor depends on the subband index since the subbands 
differ in their reaction to V g variation. The differ¬ 
ence results from the specific form of the transverse 
wave functions for separate subbands. However, we 
find that the conversion factors do not depend on the 
SO coupling constant 7 rs b- The determined values of 
the level arm for the first, second and third subbands 
C/ onv = 0.27, Cconv _= 0.26, C c 3 onv = 0.32 are used 
below for determination of the g* factors. 
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played in Fig. [4] We read out the Zeeman splitting of 
the transconductance lines in external magnetic field 
on the gate voltage AV g scale. The gl factor for k -th 
subband is then evaluated as 


1 d(AV s (B)) k 
gk fl B dB 


(7) 


The splitting of transconductance lines is plotted in 
Fig. 0 and the calculated slope determines the effec¬ 
tive Lande factors. 


Figure 3. (a) Current as given by the Landauer for¬ 
mula |5|). (b) Conductance G = dI/dV s d in units of 

Go = t-. The lines indicate the conductance steps, (c) 
The transconductance plot G = d 2 I/dV s ddV g calculated 
from (a). The green arrows indicate the voltage condi¬ 
tions for which the lowest subband follows the source or 
drain Fermi levels as the gate voltage V g is varied. The 
lines are used to determine the level arm, i.e. the gate- 
voltage to energy conversion factor, (a-c) were obtained 
for B = 0 and 7 ra b = 0. (d) The transconductance for 
B = 2.0 T. 
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Figure 4. Transconductance for x (a,d,g), y (b,e,h) and 2 
(c,f,i) orientation of the magnetic field, in the absence of 
the electron-electron interaction and SO coupling 7 rs b = 0 
(a-c), for 7 ra b = 11.4 meVnm (d-f) and 7 ra b = 40 meVnm 
(g-i). The dashed lines indicate conditions for which one 
obtains k x = 0 at the Fermi energy for a homogenous 
infinite channel with the lateral confinement adopted from 
the center of the constriction. 

As a next step in the procedure, the transconduc¬ 
tance is evaluated as a function of the magnetic field. 
The results with and without the SO coupling are dis- 


For B = (B x , 0, 0), the SO coupling leaves the split¬ 
ting of the transconductance lines AV g unchanged 
[Fig. 0a, d ,g) and Fig0a,d)[. On the other hand the 
splitting A V g for the magnetic field oriented in-plane 
but perpendicular to the current flow B = (0, B yi 0) 
is decreased by the SO coupling [Fig. 0b,e,h) and 
Fig. 0b, e)[. The effective magnetic field for electrons 
moving in the x direction polarizes their spins in the 
y direction. Therefore one should expect that SO in¬ 
teraction will leave the Zeeman effect for y orientation 
of the field unchanged. Moreover, since for B = 0 the 
spins of moving electrons are already parallel or an¬ 
tiparallel to y axis [see Fig. 0 A V g (B x ) should be ex¬ 
pected decreased by SO interaction. Nevertheless the 
calculated results are just opposite. The transconduc¬ 
tance lines for B x tends to cross at higher magnetic 
field. No anti-crossing is observed for this orienta¬ 
tion of the magnetic field, as should be expected for 
the competition between spin polarization between B x 
and the SO effective magnetic field. 

Summarizing, the results of Figure 0 In the ab¬ 
sence of the spin-orbit interaction the evaluated in¬ 
plane g* factor is isotropic for the in-plane orientation 
of the field and is nearly equal to the material con¬ 
stant applied in the Hamiltonian g = 9, and subband 
independent. The orbital effects of the perpendicular 
magnetic field B z due to the potential introduced by 
the field within the constriction enhance the evaluated 
g* factors for higher subbands also when SO coupling 
is absent. For the conversion factor which increases 
with the subband index, the orbital effects keep the 
AVg unchanged for subsequent subbands in Fig. 0c). 
Fig. § indicates, that the (i) spin-orbit interaction 
enhances the g* z factor (ii) introduces in-plane asym¬ 
metry of the g*-,g* factor (in) leaves unchanged g* 
and (iv) reduces g*. We analyze these effects in the 
next Section. 
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Figure 5. The splitting of the transconductance lines as 
read out of Fig. [4] for the magnetic field oriented in the x 
(a,d), y (b,e), and 2 (c,f) directions. The electron-electron 
interaction is neglected in these data. 



Subband index, n 


Figure 6. The calculated g* factors for various orien¬ 
tations of the magnetic field and spin-orbit coupling con¬ 
stant. The electron-electron interaction is neglected in the 
plotted results. 


IV. SO EFFECT ON OBSERVED g * FACTORS 
A. Neglected subband mixing 


In order to explain the results presented above, let 
us consider a homogenous channel with the lateral 
confinement given by a harmonic oscillator potential 

U ho = -TOeffwV- (8) 

The spatial electron wave functions in the absence of 
SO interaction for the n-th subband is of the form 


i/j n {x,y) = A n e t H n 


h 


(9) 


with normalization constant A n = i— ( m eC tu \ 4 

V 2 "n! ' irfi > 

and H n being the n-th Hermite polynomial. The 
transconductance that is used in the procedure for ex¬ 
tracting the g* factors is non-zero when one of the sub¬ 
bands enters or leaves the transport window defined 
by the source and drain potentials, hence for further 
analysis we take fc" = 0 (this condition follows the 
transconductance steps as evaluated from derivatives 
of the current - see Fig. [4]) 

We consider the effects of SO interaction for weak 
magnetic fields. We neglect the term that is quadratic 
in B z for the magnetic field oriented perpendicular to 
the plane of confinement. We adopt the spin basis 
of a- operator. The position of the bottom of sub¬ 
sequent subbands for k x = 0 is E n = Hu> (n + \). 
For a narrow channel when the value of Hoj is high 
the energy spacings between the subbands are large. 
Therefore, it seems justified to neglect the subband 
mixing m■ We diagonalize the SO Hamiltonian in 
the basis of the spin-up and spin-down basis of the 
lowest subband (n = 1 ), 

'h \ = V' ~ H a f j J ) , (10) 

For analysis we consider the Hamiltonian H (JIJ (with 
U ho term © replacing V ext ), in which the lateral SO 
field is missing due to the vanishing of the term 9 g° xt 
and k x = 0. For k x = 0, the Rashba SO interaction 
reduces to if rashba = 7rsb {cr x k y - <jy (-eB z y)}. The 
Hamiltonian matrix for the above basis 

H a , a = (<(v \H\ <t >„), (12) 


takes the form 


H 


n 


tp gn&Bz g[i&(iBy-\-B x ) 

■E'n i 2 2 

g^B(B x — iB y ) rp _ gg-^Bz 
2 2 


(13) 


From the form of the Hamiltonian it is evident, that 
for neglected subband mixing the effect of the Rashba 
interaction does not affect the Hamiltonian matrix. 
Thus the eigenvalues of the Hamiltonian above are 
same for three perpendicular orientations of the mag¬ 
netic field: 


K {(B, 0,0)) = El ((0, B, 0)) = El ((0,0, B)) 

TP L dk’BB 

= E n + a -. 


2 
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B. Subband mixing effect 


In order to describe the SO effects on the energies of 
the bottom of the subbands for the in-plane magnetic 
field orientation one needs to account for the subband 
mixing. We add the second subband to the basis 


_ "‘■ett TT 

01,^ = -'T-lC 2h H\ 

2 

_ rn eif uy 

(f> i,f = A\e 2h H\ 


m eS uj 


h 


m e ffU> 


h 


-y 


(14) 

(15) 


The Hamiltonian matrix (121 in the basis of four states 
with n = 0 and n = 1 [Eq. (101, (111, (141, and ( |l5| )] 
takes the form 



B [T] 


H o Hoi 

Hi Hi 


where H o and Hi are given by (13) and 

0 , ,rab = (m e ffuj — eh B z ) 


H oi s 


. y/2hm ett ui 


^/2 hm, 

(m e jjoj + ehB z ) 0 


Note, that for B z = 0 the SO constant 7 rs b gets ef¬ 
fectively enhanced for a larger values of Hco. Now, for 
B = ( B Xl 0, 0) the eigenvalues are 


E n,a = KU - 


(- 1 )” 2m eS ^l h 


/t \2 

(M +a —-— 


The Zeeman splitting for a given n is still not affected 
by the SO interaction - in agreement with Fig. [6] For 
B = (0, By , 0) one finds 


E n,a = ^ - 


(- 1 )" . / 2 TO cff a; 7 r 2 sb 

2 V h 


(hu + aguBByY. 

(16) 


For nonzero 7 rs b the spin splitting is reduced [see Fig. 
[7], in agreement with the results of Fig. [g] For B = 
(0,0 ,B Z ) one obtains 


El „ 


(-1 r 


2Weff7rsb 

h 3 CJ 


hjJ 


(hu + agfj, B B z ) + 
aeH 2 B 7 . 


TO eff 


(17) 


For the H z magnetic field the SO interaction is more 
pronounced when the confinement strength is smaller 
(hco —> 0). 


Ene rgy lev els as calculated in the basis of Eq. 
1 14|), and (151 for hw = 2meV. Black dashed 


Figure 7. 

( 10 [), |TT}, _ _ 

lines indicate the results obtained in the absence of the SO 
coupling. The results are identical to the ones obtained for 
B = (.Eh,, 0,0) in the presence of the SO coupling - with 
a precision of a B-independent energy shift of the energy 
levels. The red lines show the results for B = (0,0, B z ) 
and yrsb = cx 3 dF z = 10 meVnm. The light and dark blue 
lines indicate the results for B = (0,B y ,0) for -y rs i, = 10 
meVnm and 7 ra b = 40 meVnm, respectively. A constant 
shift of the energy levels for different values of y rs b was 
applied. 


The numerical data are given in Fig. [7] for hui = 
2meV. The calculation reproduces all the features ob¬ 
served in Fig. [6] listed at the end of the previous sub¬ 
section. 

We conclude that, the experimental procedure for 
extraction of g* factors from transconductance de¬ 
tects the bottom of the subbands for which k x = 0 . 
For k x = 0 and B y = B z = 0, 07 commutes with 
the Rashba and total Hamiltonian. Therefore, for 
B = (B x , 0,0) the Hamiltonian eigenstates remain 
eigenstates of a x operator even in the presence of the 
SO interaction. Hence, for the magnetic field oriented 
in the B x direction the SO coupling does not influ¬ 
ence the spin states nor the g* factor. The crossing 
of transconductance lines of Fig. [ 4 ] for B x follows (see 
Fig. 4}g)). Hence, the spins at the bottom of the sub¬ 
bands get oriented by the SO interaction in the ±x 
direction and not in ±y direction as expected from 
the orientation of the effective SO magnetic field. For 
By 0 the splitting of energy levels is reduced by 
the SO interaction. Moreover, the description of the 
SO effects at the bottom of the subbands and the re- 















































suiting anisotropy of g* factor requires account taken 
for the subband mixing. The approximation of a sin¬ 
gle subband [22. is not suitable for evaluation of the 
g* factors modification by the SO interaction as cal¬ 
culated by the main simulation of this paper, and as 
observed in the transport experiments [23] EU [33] [35]. 
Moreover, for Ep between and E 2 the SO coupling 
effects involve coupling of a transmitting (n = 1) and 
evanescent (n = 2) modes, i.e. the subband mixing 
is present also when one of the subbands is above the 
Fermi level. 


V. INTERACTION EFFECTS 

For evaluation of the interactions effects on the 
g* factors we used the density function theory [45] . 
The calculations for the interacting systems were per¬ 
formed for somewhat smaller computational box of 
length L x = 600nm and width L y = 160nm with grid 
spacing Ax = 4nm, and N = 150 electrons inside the 
2DEG and the QPC confining potential given by a 
similar formula as in Eq. ([2]) 

Kxt = V gllte (x, y\ 300nm, 0, 50nm, 34nm) 

+ Vg.nte(x, y ; 300nm, 160nm, 50nm, 34nm). 

The approach used in this paper is a development of 
the one used in Ref. m , with an extension to spin- 
orbit coupling effects. The details are given in Ap¬ 
pendix B. The DFT is used to evaluate the Fermi 
energy, and the spin-dependent potential landscape. 
The calculation is performed as a function of V a and 
B for V s d — 0. For larger V s d we introduce the bias as 
an additional linear term as in Section IIIIBI 

The transconductance plot with the electron- 
electron interaction included is given in Fig. [8] for 
B = 0. The electron-electron interaction introduces 
asymmetry in the transconductance lines as functions 
of the gate-voltage, clearly visible for the first conduc¬ 
tance step. A different path of the lines for V g below 
or above the QPC conductance step is related to the 
appearance or disappearance of the electron density 
related to a specific subband depending on the gate 
voltage [Fig. §a)[. For V g > 250 meV the electron 
density is removed from the center of the QPC [Fig. 
[9jc)[ and the screening of the gate potential by the 
electron gas is missing, hence a larger slope of po¬ 
tential in Fig. §b). Similar effects of the pinch-off 
for the screening were obtained in Ref. EZI, in which 
the transconductance was not disussed. The present 


paper indicates their consequences for asymmetry of 
the transconductance lines. The asymmetry of the 
transconductance is the strongest for the first sub¬ 
band [Fig. [8], but appears also for higher subbands, 
which is accompanied by deviation of the potential 
landscape from a linear dependence on V g [Fig. |9][b)]. 
In the experiment [23j the transconductance lines are 
(i) asymmetric with V g , (ii) get steeper at the higher 
absolute value of the gate voltage side, and (in) the 
asymmetry is reduced for larger n. Our calculation 
reproduces all these features and explains that their 
origin is the screening of the gate potential by the 
two-dimensional electron gas. The evaluated conver¬ 
sion factors are C( onv = {0.14,0.10,0.12} for the first, 
second, and third subbands, respectively. Note, that 
the reduction of the conversion factors and their vari¬ 
ation with the subband index is the effect of the gate 
potential screening. 



150 200 250 

V„ [mV] 


Figure 8. Transconductance for 7 rs b = 0 as obtained with 
the electron-electron interactions described by DFT. 


The transconductance dependence on the magnetic 
field results is given in Fig. [TO] with the positions 
of the maxima marked by blue shadows. The evalu¬ 
ated effective g* factors are displayed in Fig. 11. The 
features introduced by the electron-electron interac¬ 
tion are revealed by comparison with Fig. 6 of Sec. 
IIIB We find (i) An enhancement of the g* above 


the nominal value g = 9 for the first n = 1 subband 
index for each direction of the magnetic field, (ii) 
The largest enhancement (a factor of 3) is found for 
perpendicular magnetic field in case B = (0,0, B z ). 
In the absence of SO interaction the in-plane Lande 
factors g* and g *, exhibit the same behavior as a func¬ 
tion of n. When the Rashba SO in introduced with 
7 rs b = 11.4meVnm the value of the g( drops faster 
with n than g* leading to an anisotropy for in-plane 
magnetic field. The results, including the last feature 
remain in a good agreement with the experiment of 
Ref. [23 . The g* / g* anisotropy is weak in Ref. [23]. 
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V g [mV] 



100 150 200 25050 100 150 200 250 J '%0 
V [mV] K [mV] 


100 150 200 250 
V [mV] 


Figure 10. The transconductance as a function of mag¬ 
netic field, for (a) B = (B x ,0,0) (b) B = (0,B y ,0) and 
(c) B = (0,0, B z ) for 7 ra b = 11.4nm 2 . The insets show 
schematically the orientation of magnetic field for each 
case. The blue lines indicate the positions of the transcon¬ 
ductance peaks used for the evaluation of the g*. 


Figure 9. (a) The conductance as a function of the gate 
voltage for V a d = 0 corresponding to[8|a). The effective 
potential obtained from self-consistent DFT calculations 
(b) and the electron density (c) at the center of QPC con¬ 
striction. 



Subband index, n 


A stronger anisotropy should be expected to appear 
in InSb quantum point contacts m- 


For larger value of 7 rs b = 40meVnm (i) the effec¬ 
tive g* is almost the same as for other Rashba pa¬ 
rameters, with a small but visible shift towards larger 
values, (ii) g* drops in a more distinct manner with 
n. This in an effect encountered already in the pre¬ 
vious subsection for the model without the electron- 
electron interaction. For the out-of-plane magnetic 
field we observe a strong enhancement of the g* fac¬ 
tor with respect to the in-plane g* and g* factors in 
agreement with Ref. [25] ■ We obtain more or less a 
constant behavior in n for the 7 rs b = 0. For larger 
SO coupling constant the g* tends to grow with n. In 
the experiment [25] the g* drops with n. In Fig. H c ) 
we have shown that even without electron-electron in¬ 
teraction the orbital effects may lead to increment of 
g* for higher subbands. Additionally, those effects de¬ 
pend highly on the geometry of the sample thus, one 
may expect, that the effect will differ from sample 
to sample. On the other hand the results for g*, g* 
are free of orbital effects and the results should not 
be sample dependent, hence good agreement of the 
present results to the experimental data. 


Figure 11. The calculated values of \g*\ for three different 
orientations of the magnetic field (a-c) and for different 
strengths of the 7 ra b as a function of the subband index. 
The error bars denote the standard deviation of the linear 
fit performed for AV g (B) in Fig. 10 The dashed line on 
each plot indicates the nominal value of the g = 9 used 
in simulations. In case of B = (0, 0, B z ) we used data for 
the range of the magnetic field for which A V g (B) had the 
linear behavior. 


VI. SUMMARY AND CONCLUSIONS 

We described a numerical model to solve the trans¬ 
port problem for a biased InGaAs quantum point con¬ 
tact for a two-dimensional electron gas with a confine¬ 
ment potential inducing narrowing of the conducting 
channel. We studied the effects of the Rashba spin- 
orbit and electron-electron interaction for transcon¬ 
ductance in the regime of fractional conductance 
quantization in external magnetic field. The calcu¬ 
lated transconductance was used to extract the effec¬ 
tive Lande factors as functions of the external mag¬ 
netic field orientation. For the extraction we simu¬ 
lated the standard experimental procedure using the 
gate-voltage to energy conversion factor evaluated 
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from the transconductance and its reaction to the 
source-drain bias. 

We discussed the anisotropy of the g* factors look¬ 
ing for the effects due to the spin-orbit interaction 
and the interaction effects. We have established that 
the spin-orbit interaction alone enhances the g* factor 
for the out-of-plane orientation of the magnetic field. 
Moreover, the Rashba interaction introduces the in¬ 
plane anisotropy of the electron g* factors. The g* fac¬ 
tor for the magnetic field parallel to the current flow 
direction is unaffected by the SO interaction, which 
reduces the g* factor for the perpendicular orienta¬ 
tion of the in-plane magnetic field. We explained that 
this effect - counterintuitive from the point of view 
of the effective SO-related effective magnetic field 
is due to the subband mixing for the bottom of the 
subbands entering the transport window, for which 
non-zero transconductance is found. The explanation 
was based on an analytical perturbation analysis. 

The electron-electron interaction alone leads to a 
pronounced enhancement of the g* factor for out-of- 
plane magnetic field. Moreover, the screening of the 
gate potential by the electron gas as described by the 
DFT approach induces asymmetry of the transcon¬ 
ductance lines on the gate potential - bias plane that 
is most pronounced near the pinch-off of the trans¬ 
port across the quantum point contact. For combined 
electron-electron and Rashba interaction effects the 
calculated g* factors agree qualitatively with the ex¬ 
perimental data. We found that the electron-electron 
interaction preserves a weak dependence of the g* fac¬ 
tor on the subbands index for the magnetic field ori¬ 
entation parallel to the current flow, and the drop of 
the g* factor for the other perpendicular in-plane ori¬ 
entation of the magnetic field vector. 
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APPENDIX 

A. Description of the scattering problem 

In order to solve the scattering problem of the Fermi 
level electrons inside the device discussed in this paper 
we implemented the QTBM method in the finite dif¬ 


ference formalism formally similar to the tight-binding 
(TB) approach. The procedure developed here is a 
generalization of our previous method [IS] with in¬ 
clusion of the spin-orbit coupling. We work with 
the wave-function description of the transport with 


spinors component x = 


coupled by the SO in¬ 


teraction. The derivation is similar to the wave func¬ 
tion matching method |19| . The general Schrodinger 
equation for the system presented in Fig. [12] can be 
written in a compact form 


'Tn—lXn—l T H n Xn T ' J ~nXn+l — Xn-> (18) 


where T n ± i is the coupling matrix between slice n and 
nil and H n is the TB Hamiltonian describing n-th 
slice (see Fig. 121. The total wave function in the 


device contains contribution for all the slices 


(u,v) = y iX k(v)S u , k . 


(19) 


k=1 



Figure 12. (a) Schematic representation of the tight- 

binging device considered in this paper. The system can 
be viewed as ID chain consistent of slices \k coupled by 
the matrix r. Here y is represented by the vector of size 
equal to number of points in transverse direction. The 
left and right leads are assumed to be ideal semi-infinite 
contacts. In our case points on the mesh are separated by 
a distance A*. The whole system consists of N x N y grid 
points, (b) Coordinates of the scattering wave function on 
a square TB mesh. 
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In order to find the scattering wave function one 
has to apply proper boundary conditions at the de¬ 
vice/lead interface which will satisfy the continuity of 
the wave function and its derivative at the interface. 
The standard approach for this problem is to find a 
general solution in the homogeneous semi-infinite lead 
(r„ = t 0 for any n in the lead - see Fig. 12). Thus 


the Eq. (18 1 for the left lead can be written as 


ToXn-l + HoXn + T' 0 Xn+l = E F Xn- (20) 


This equation can be solved by applying Bloch substi¬ 
tution Xn = A n x which assumes that the wave func¬ 
tion between two nearest slices differs by a constant 
factor A. The number A can be interpreted as phase 
which the wave function gains along a distance Ax 
between two slices and can be written in a plane wave 
form A = e lkAx . After this substitution and division 
by A n_1 the Eq. ( |20| ) takes the form 

ToX + H 0 Xx + TqA 2 x = E f Xx, 


which with substitution x = Xx can be transformed 
into the generalized eigenvalue problem for {x, x} and 
eigenvalue A 


0 I 

T o Ho ~ 




( 21 ) 

Upon solution of this eigenequation one obtains the 
transverse modes Xm in the lead and the correspond¬ 
ing wave vectors k m . We solve the eigenproblem (211 
numerically with the LAPACK library [SUJ ■ In order 
to distinguish between the traveling and evanescent 
mode we use fact that for the former |A m | = 1 and for 
the later: |A m | > 1 is decaying reflected mode present 
in the left lead near the scatterer - and |A m | < 1 is the 
decaying transmission - present in the right lead close 
to the scatterer. We distinguish the traveling modes 
Xk between the incoming {Afc j+ , Xk,+ } an d outgoing 
{A k,~ , Xk ,-} ones modes by the sign of the carried cur¬ 
rent flux <t>k’- for (f>k > 0 the mode was identified with 
Xk. + and when (f>k < 0 => Xk,-- For the evanescent 
modes the flux is always zero. The formula for </>& can 
be obtained by considering the continuity equation 

dt 3 ’ 

in the TB formalism which leads to matrix formula 
for the current flux 0k oc 3 f At-v! rtyt-J [51]. We de¬ 
note M as the number of the traveling modes in the 


lead. Note that the number of the evanescent modes 
in general is infinite but when we deal with TB like 
Hamiltonian this number is always finite and equal to 
M e = N — M, where N is the size of the vector x- 
Having sorted transverse modes one may write gen¬ 
eral solution for the left semi-infinite lead at position 
(u, v ) on the mesh 


M M+Me 

Vl(u,v) = ^a fei+ A£ )+ Xfc,+ + r k- X k,-Xk-, 

k k 

( 22 ) 

where we assume that AJ1 _Xk,~ for k> M + 1 corre¬ 
spond to decaying evanescent modes and ak,+ is the 
scattering amplitude of the incoming mode and Tk,~ 
is the reflection amplitude. Analogical equation can 
be written for the right lead 

M M+M e 

*n(u,v) = Y,ak,- x l-Xk,-+ E r k,+ X k,+Xk,+ - 
k k 

We will denote by Xfc.± i v ) the v component of the 
Xk,± vector. The idea of the QTBM is to match the 
wave functions inside the leads \l/n and 'Fl with the 
\Fd. The continuity condition at the left lead interface 
implies 


'F l (w = 0 , v) = Td(u = 0 , v) = xo, = ^d(w, v). 


Thus from Eq. (22) we have 

M 


M+M e 


^d(w = 0,v) = Xo(v) = J^a ki+ Xk,++ E r k.-Xk,-- 


k k 

Following refs [42, 43), 48; [52] we express rk,~ in terms 
of 'F l (0,v) and cik,+ 


m+m b 

r k,-= E S kp-(Xp,-\Xo) 

p =1 

M+M e M 

- E E &kp,— -Apq,— Clq,-\- 5 (23) 

p=l q= 1 

with (Xp-lxo) = AzXEiXp-MxoCf) being the 
inner product, = (Xfc,-|Xp,-) and A k P - = 

(Xfc,-IXp,+) are the transverse modes overlap matri¬ 
ces. Note that S- is always a square matrix of di¬ 
mensions (M + AI e ,M + M e ), while has dimen¬ 
sions (A/ + M e , M) and is a square matrix only when 
one neglects the evanescent modes in the calculations. 
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The second condition that wave function at the in¬ 
terface lead/device should fulfill is the continuity of its 
derivative at this point. Numerically it can be written 
in the following way 


'M+l, v) - v) = >F +M - (24) 

= Xi(u)-X-i(u), (25) 

which is the finite difference equivalence for the ana¬ 
lytical condition 


0 %, 

dx 

The left side of the 


Eq. (221 


x—0 

Eq. 




dx 


x—0 


(24) can be evaluated using 


M 

T L (+l,u) - T l (-1,u) ='^ak,+A ki+ Xk,+ 

k 

M + M e 

+ Y r k-Ak-Xk- 

k 

=Xi{v) ~ X-i(u) 


where A kyd = Xk,d - and d = (+,—). Now in¬ 
serting the expression for r k ,- from Eq. ( |23| ) to the 
equation above leads to 


N 

Xi(v) - X-i{v) = T+ + A~iXo(i), (26) 

i —1 

with 

M M+M e 

='Yh a k,+A k ,+Xk,+{v) ~ Y & k,- A k,-Xk,-(v) 

k k 

(27) 

M+Me M 

O^k,— = ^ ^ ^ ^ &kp,— -Apq : — Q>q,-\- (28) 

P= 1 <7=1 

M+M e M+M e 

A v,i = Ax Y Y A k-Xk-{v)S kp -xl-(i)- 

fc=l p=i 

(29) 


With Eq. (261 we can calculate the value of the wave 
function outside the device X-i- Using the vector no¬ 
tation (see Fig. 12 and Eq. @) we get following 
formula 


X-i = Xi - r+ - A xo, 


which we put to the Eq. ( fl~8j ) for n = 0 in order 
to remove the reference to the slice %-i outside the 
device 


T— 1 (xi - T + - A xo) + H oxo + T oXi = EfXo, 


where we assume, for simplicity that lead are far from 
the scattering center, which is usually located in the 
center of the device, thus we can put t_i = To. Sim¬ 
plifying the equation above we get expression for the 
boundary condition in the left lead 


{H 0 - IE f - r 0 A ) xo + (to + rj) Xi = t 0 F + . 


Analogically to Eq. (261 we get similar equation for 
the right lead boundary condition 


N 




N+l,v 


— ’Fn-i,d = r 




where the definitions for vector T. u and matrix Ab 
can be obtained from equations ( 27][29 1 by changing 
sign H—> — and-> +, and noting that the calcula¬ 

tions of those quantities are performed for A-th slice 
of the device. From the equation above we calculate 
the formula for N + 1 slice using the vector notation 


xn+i = xn-i + r + a + xn, 


and we substitute it to the Eq. (181 for n = N assum¬ 
ing that tn = Tn-Ij 


tn-iXn-i+-HnXn+t^ (xn-i + T + A + xn) = -EfXn, 
which leads to final formula for the right lead 

(r N + tJj) Xn-i+(-Hn - IE f + t^A + ^ xn = -r^r - . 

Finally the whole matrix equation can be written in 
the form of the system of linear equations A\ = b, 
with block matrices 

( H' 0 -t 0 A~ t 0 + tJ ■ 0 \ 

To H[ t{ 


A = 


V 


0 


r t H' 

r N —2 “ ™ 


r t 

' N— 1 


tn+tI H' n +tLa + ) 


(30) 


with H’ n = H n -IE f , X = (Xo, Xii • • • > Xn-i, Xn) 
and b = (toF + , 0,..., 0, —rj^T - ). In order to obtain 
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the scattering wave function for the p-th mode incom¬ 
ing to device from the left one must put (ik.+ = <5fc,p 
in the definition of the T + and af. = 0 for all k in 
the r-. For the current flow from right to left we put 
ak t - = Sk tP in r- and <ik,+ = 0 in the definition of the 
r+ . We use the PARDISO library to solve the system 
of linear equations [531 . 

Let us denote the k-th scattering mode by where 
k = (1,2enumerates modes in the left lead. 
Once the system of linear equations is solved one may 
calculate the reflection amplitudes from Eq. (231 by 


calculating the value of ru = and transmission 

amplitudes tfc = Tk,+, which in case of the transverse 
mode incoming from the left lead we have 


M+M e 

tk = ^fc,+ = 'y ' Skp,+ (Xpj+I^N,!;) • 

p= 1 

Using calculated values of rk and tj~ we calculate the 
reflection R p and transmission T p probability of the 
p-th incoming mode from the left lead with equation 


M 


k =1 
M 


k —1 


</>L ,k- 


0L,p,-f 

^R,fc,+ 




where 4>{l,p}. fc,{+,_} denotes the quantum flux which 
mode number k is carrying in the L - left lead or R - 
right lead in the +- right direction or - left direction 
(outgoing modes in the left lead). From that we calcu¬ 
late the conductance of the system using the Landauer 
formula 



p=i 


The tight-binding form of the Schrodinger equa¬ 
tion was obtained by applying the gauge-invariant 
kinetic-energy finite difference discretization |54| of 
the Hamiltonian given by Eq. (JT|). This discretiza¬ 
tion leads to following matrix equation for the device 
wave function 

( yf T A 1 " yf si + \ ( yf \ 

v xL s I + r lv + A iv ) F v xL J ’ 

(31) 


with 


D a = 
uv Ax 2 

rj-ida _ 

^ uv 


2 1 

+ eVext (u, v) + -agn B B z 


1 


S a = ~d 


2Ax 2 

^Orsb 


D +ido-7iat-Ex(^,'y)Ax 


= -ds 


2 Ax ’ 

*7rsb -diAx 2 B z 


2Ax 


S b — A 1 

A CT = \- a T~ a A- v CT T +a A- T~ a 

uv Auv— 1 uv— 1 1 Auv+1 uv 1 Au—lv'u— Id 

, V <T +(T 

r Au+lv ' uv ? 

I"uv = Xu-lv'S'v a + Xu+lv^ CT + Xuv-lS 


+ Xuv+1^ + ! 


T da = - 


2Ax 2 


0 —diAx 2 B z v—idaj\ a tB y (u,v)Ax 


where d = ±, and the electric fields generated 
by the in-plane external potentials E x (u, v) = 
(V ext {u+l 1 v)-V ext (u-l 1 v))/2Ax and E y (u,v) = 
(Text (u,V+ 1) - Uext (u,V~ 1)) /2Ax. 

We have mapped the 2D+spin coordinates Xuv into 
ID coordinate using following formula: w{u,v,a) = 
( u ^ T y + v) hr order to get row or column index 
of the TB matrix. Note that such parametrization of 
w divides the entire matrix into two major blocks 


/ H t S n 

V S' 4 - 1 H l 


e f 


X 1 

X 4 


with H a being the Hamiltonian of the upper or lower 
part of the spinor and S r7rJ being the coupling ma¬ 
trices. When there is no SO interaction or magnetic 
field is oriented only in z direction the matrix can be 
easily divided into two independent subspaces, twice 
times smaller. In order to calculate the conductance 
we applied the QTBM described above to Hamilto¬ 
nian Eq. (311. The dispersion relations presented in 
Fig. i were obtained by substituting the plane wave 
of form 


_ ikuAx I Xl 

Xuv - e l I 

Ad 


into Eq. (311 and then calculating the eigenvalues of 


created eigenproblem as a function of k. 
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B. Description of the many electron problem have calculated the electron densities 


We solve the many electron problem using the 
Density Functional Theory (DFT) with Local Spin- 
Density Approximation (LSDA) for the evaluation of 
the exchange-correlation functional. The approach 
used in this paper is similar to those described in our 
previous work m or in [56]. All the DFT calculations 
were done for Vsd = 0 case. The Kohn-Sham (KS) or¬ 
bitals were obtained by diagonalizing the Hamiltonian 
given by Eq. (311 with the spin dependent exchange- 
correlation potential K| c and Hartree Vr potential 
added to the diagonal terms 


D 


S 

UV 


= ^2 + eKxt {u, v ) + -sgfi B B z 
+ eV H (u,v) + eFxc(u,D). 


(32) 


The periodic boundary conditions were applied in the 
x direction in order to simulate open infinite leads. We 
use the Attaccalite parametrization for the !/£<-, poten¬ 
tial 156], which depends on the spin-up electron den¬ 
sity n- 1 - and spin-down n±. The Hartree potential was 
calculated by numerical integration of the Coulomb 
interaction term 


g 

Vh(x, y) = —— (-J d (x, y\ 20nm) + I e {x, y; 0)), 
47T€€o 


where we used e = 12.4 for InGaAs, and Id is the 
coulomb integral over the positively charged donor 
layer and I e - integral over the total electron density 
in 2DEG layer and 


I x {x,y,z x ) 


dx'dy' 


n x {x\y') 

\(x,y,0) - (x 1 , y', z x ) \' 


In order to simulate the infinite leads we performed 
the integration including 8 copies of the system in 
left and right direction. In the formula above we as¬ 
sume that the 2D donor layer is located = 20nm 
above the layer of 2DEG (z e = Onm) and the den¬ 
sity of the ionized donors in that layer is constant. 
The total number of donors was = 150 in our 
calculations and the total electron density was cal¬ 
culated as n e (x,y) = n^(x,y) + ni(x,y). In the 
first step of DFT iterative scheme we assumed that 
n e = n-j- = n± = 0, then we performed the numer¬ 
ical diagonalization of Hamiltonian ( 31|32 ) with the 
FEAST eigenvalue solver ML calculating the first 
lowest M e = 500 eigen-energies ef and correspond¬ 
ing KS spin-orbitals xf — Xi( u ,v). From that we 


M e 


= n^u 1 v) = Y J f{^EY,T)\xl\ 2 , (33) 


k =1 


where / is the Fermi-Dirac distribution 


f(4-,E Fi T) = 


1 




and the Fermi energy Ep was obtained from the charge 
neutrality condition 


M e 


^d = E/(4;^F,T) + /(4;£; F ,T). (34) 


k =1 


We solved both Eq. (331 and (34) for densities n a 
and Ep using the bisection method. Once the n e ,n ^ 
and nj. were obtained we have calculated the Vr and 
K| c potentials then we solved the KS Eq. ( 31|32 > 
again. The whole procedure was repeated until the 
density stopped to change significantly between two 
consecutive iterations. 

In order to improve the convergence of the DFT 
iterative scheme we have implemented a) the temper¬ 
ature annealing b) density mixing method. 

In case of temperature annealing we start the cal¬ 
culation with temperature T =10K, then after the 
convergence was obtained we use calculated electron 
densities as input densities for lower temperature cal¬ 
culations. The annealing was performed in three steps 
for T = {10AT, 4AT, 0.5AT} temperatures. 

In case of density mixing we applied the simple lo¬ 
cal and linear mixing scheme [58j. Let us denote by 
4 the electron density with spin er in k- th DFT itera¬ 
tion, and by n„ +1 = F(n^) electron density obtained 


after diagonalization process from Eq. (331. If we 
define the residual vector R(n £) = F(n ^) — n* then 
standard approach for almost all mixing schemes start 
from equation 


4+ 1 := n\ + (35) 

where JA = i s the Jacobian matrix, Ri(n Jj) 

is the i-th component of the residual vector (the 
residuum value for i-th point in space) and := repre¬ 
sents the numerical substitution. In the simplest case 
one may assume that Jacobian is a constant diagonal 
matrix J a = — al , where 0 < a < 1 is the mixing pa¬ 
rameter, chosen arbitrarily for a given problem, then 
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the mixing formula (351 for electron density in the 
next steps takes the from 


n k+1 :=n k +a (F{n k ) - n k ) = an k+1 + (1 - a) n k . 

(36) 

In the second approximation one may assume that 
Jacobian is a diagonal matrix, which change during 
the iteration, of form 

[JkV 1 = -diag«,a^,...,a^), (37) 


af = a max ) otherwise set af = a m i n . Note that the 
method works for k > 3, hence for k = 1 we use the 
direct substitution n k+1 := F(n k ). In second itera¬ 
tion k = 2 we used the simple mixing formula (361 
with a = 0.05. Then for k > 3 we implemented the 
local mixing algorithm given by Jacobian (371, with 
a m ; n = 0.05 and a max = 0.8. According to the pa¬ 
per |58| this method is almost as efficient as the more 
complicated Broyden method 1521 - 


where N is the size of the vector R. This method 
start with af = af = ... = a a N = a m ; n then at each 
iteration the parameters af are changed according the 
following algorithm: if J?j(n^)iZj(n* _1 ) > 0 then in¬ 
crease the value of af by a m i n (if af > a max then 


Once the convergence of the KS equations was ob¬ 
tained we use the obtained from Eq. (331 and (34 ) 
Fermi energy Ep and calculated Vh and I4|c poten¬ 
tials in the one electron picture method described in 
section IVI Al 
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